knitr::opts_chunk$set(echo = TRUE)
#install.packages("pacman")
pacman::p_load(here, dplyr, ggplot2, ggpubr, viridis, ggh4x, stringr,
forcats, tidyr, lme4, emmeans, plotly)
# Ryssa's go-to theme
theme_set(theme_light()+
theme(
plot.title = element_text(size=rel(1.2),face="bold"),
axis.title = element_text(size=rel(1),face="bold"),
axis.text = element_text(size=rel(1),colour = 'black'),
strip.text = element_text(size=rel(1),colour = 'black',
face = "bold"),
legend.text = element_text(size=rel(1)),
legend.title = element_text(size=rel(1),face="bold"),
panel.grid = element_blank()))
Pairs of people played the mirror game. We measured how well people were able to mirror each other’s spontaneous movements to understand how movement synchrony influences response inhibition and the related patterns of brain activity.
For this example, we keep the columns:
ID participant identifierleader whether the participant (P) or helper (C) was
leading the gamemeanSim the mean level of movement synchrony during the
gamegroup whether the participant belonged to the
synchronised or control groupposeSimilarity <- read.csv(here::here("data/pose_similarity.csv")) %>%
filter(ID <500) %>%
select(-c(cutoff, nFrameskip, nFrames, sdSim)) %>%
mutate(group = factor(case_when(ID <200 ~ "Synchronised", TRUE ~ "Control")),
leader = factor(case_when(leader == "C" ~ "Helper",
leader == "P" ~ "Participant")))
head(poseSimilarity)
## ID leader meanSim group
## 1 101 Helper 0.7625684 Synchronised
## 2 101 Participant 0.8129776 Synchronised
## 3 102 Helper 0.7452378 Synchronised
## 4 102 Participant 0.8203161 Synchronised
## 5 103 Helper 0.8224377 Synchronised
## 6 103 Participant 0.7534514 Synchronised
While people played the mirror game, we tracked the level of synchrony and the complexity of the mirror game movements. The participants also filled in questionnaires about themselves. We were particularly interested in embodied movement knowledge. One such measure is “body competence”, which indexes how strongly people believe that their body can accomplish physical activities.
Moffat et al. (2024). Dyadic body competence predicts movement synchrony during the mirror game.
For this example, we keep the columns:
ID participant identifierM_similarity mean movement synchronyM_entropy mean movement complexity computed using
entropymean_BCQ mean body competence score per dyadRole whether the participant was leading the game
(leader) or following the helper (following)We also create a new column:
BCQ_half to categorise dyads as being above or below
the median of mean_BCQposes_traits <- read.csv(here("data/mirrorring_embodiment.csv")) %>%
select(ID, M_similarity, M_entropy, BCQ, Conf_BCQ, Role, video_code)
# Calculate mean and z-scores for helper and participant BCQ
poses_traits1 <- poses_traits %>%
# Calculate the mean separately for each row (pair)
rowwise() %>%
mutate(mean_BCQ = mean(c(BCQ, Conf_BCQ)), # BCQ = body competence
BCQ_half = factor(case_when(mean_BCQ > 7.5 ~ "High",
TRUE ~ "Low"))) %>%
select(ID, M_similarity, M_entropy, mean_BCQ, BCQ_half, Role)
head(poses_traits1)
## # A tibble: 6 × 6
## # Rowwise:
## ID M_similarity M_entropy mean_BCQ BCQ_half Role
## <int> <dbl> <dbl> <dbl> <fct> <chr>
## 1 114 0.888 0.0555 9.5 High Following
## 2 114 0.907 0.149 9.5 High Following
## 3 114 0.629 0.107 9.5 High Following
## 4 114 0.841 0.0384 9.5 High Following
## 5 114 0.735 0.0637 9.5 High Following
## 6 114 0.557 0.0874 9.5 High Following
Let’s compare the level of synchrony observed during the mirror game
(group = Synchronised) vs. the control game
(group = Control), and find out if it matters who was
leading (Leader or Helper in leader).
head(poseSimilarity)
## ID leader meanSim group
## 1 101 Helper 0.7625684 Synchronised
## 2 101 Participant 0.8129776 Synchronised
## 3 102 Helper 0.7452378 Synchronised
## 4 102 Participant 0.8203161 Synchronised
## 5 103 Helper 0.8224377 Synchronised
## 6 103 Participant 0.7534514 Synchronised
model <- lmer(meanSim ~ group*leader + (1|ID), data = poseSimilarity)
summary(model)
## Linear mixed model fit by REML ['lmerMod']
## Formula: meanSim ~ group * leader + (1 | ID)
## Data: poseSimilarity
##
## REML criterion at convergence: -219
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.46578 -0.39503 -0.02317 0.46849 2.29187
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.002627 0.05125
## Residual 0.005684 0.07539
## Number of obs: 120, groups: ID, 60
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.46831 0.01664 28.137
## groupSynchronised 0.31627 0.02354 13.436
## leaderParticipant -0.07133 0.01947 -3.664
## groupSynchronised:leaderParticipant 0.11021 0.02753 4.003
##
## Correlation of Fixed Effects:
## (Intr) grpSyn ldrPrt
## grpSynchrns -0.707
## ledrPrtcpnt -0.585 0.413
## grpSynchr:P 0.413 -0.585 -0.707
model %>% emmeans(specs = ~ group + leader)
## group leader emmean SE df lower.CL upper.CL
## Control Helper 0.468 0.0166 105 0.435 0.501
## Synchronised Helper 0.785 0.0166 105 0.752 0.818
## Control Participant 0.397 0.0166 105 0.364 0.430
## Synchronised Participant 0.823 0.0166 105 0.790 0.856
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
model_estimates <- model %>% emmeans(specs = ~ group + leader) %>%
as.data.frame()
head(model_estimates)
## group leader emmean SE df lower.CL upper.CL
## Control Helper 0.4683098 0.01664421 105.46 0.4353091 0.5013105
## Synchronised Helper 0.7845749 0.01664421 105.46 0.7515742 0.8175756
## Control Participant 0.3969787 0.01664421 105.46 0.3639780 0.4299794
## Synchronised Participant 0.8234531 0.01664421 105.46 0.7904524 0.8564538
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
model_estimates %>%
ggplot(aes(x = group, y = emmean, color = leader))+
geom_point(size = 4, position = position_dodge(0.7))+
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(0.7), width = 0.1)+
scale_color_manual(values = c("#89398a", "#398a89"), name = "Leader")+
ylim(0,1) +
labs(x = "Group", y = "Parameter estimate and 95% CI")
### Geom_bar+point()
estimates <- model_estimates %>%
ggplot(aes(x = group, y = emmean, fill = leader))+
geom_bar(stat = "identity", alpha = 0.9, position = position_dodge(.9))+
geom_point(size = 4, position = position_dodge(.9))+
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(.9), width = 0.1)+
scale_fill_manual(values = c("#89398a", "#398a89"), name = "Leader")+
ylim(0,1) +
labs(x = "Group", y = "Parameter estimate and 95% CI")+
theme(legend.position = c(0.25, 0.85))
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
estimates
model %>% emmeans(specs = ~ group) %>% pairs()
## NOTE: Results may be misleading due to involvement in interactions
## contrast estimate SE df t.ratio p.value
## Control - Synchronised -0.371 0.0191 58 -19.449 <.0001
##
## Results are averaged over the levels of: leader
## Degrees-of-freedom method: kenward-roger
group_contrast <- model %>% emmeans(specs = ~ group) %>%
pairs() %>%
as.data.frame()
## NOTE: Results may be misleading due to involvement in interactions
model %>% emmeans(specs = ~ leader) %>% pairs()
## NOTE: Results may be misleading due to involvement in interactions
## contrast estimate SE df t.ratio p.value
## Helper - Participant 0.0162 0.0138 58 1.179 0.2433
##
## Results are averaged over the levels of: group
## Degrees-of-freedom method: kenward-roger
leader_contrast <- model %>% emmeans(specs = ~ leader) %>%
pairs() %>%
as.data.frame()
## NOTE: Results may be misleading due to involvement in interactions
contrast_df <- rbind(group_contrast, leader_contrast) %>%
mutate(contrast = case_when(contrast == "Helper - Participant" ~ "Helper\n-Participant",
contrast == "Control - Synchronised" ~ "Control\n-Synchronised"))
head(contrast_df)
## contrast estimate SE df t.ratio p.value
## 1 Control\n-Synchronised -0.37136976 0.0190943 58 -19.449245 4.311649e-27
## 2 Helper\n-Participant 0.01622648 0.0137647 58 1.178847 2.432725e-01
contrasts <- contrast_df %>%
ggplot(aes(x = contrast, y = estimate, fill = contrast))+
geom_hline(yintercept = 0, color = "grey", linetype = "dashed")+
geom_errorbar(aes(ymin = estimate-(SE*1.96), #*1.96 for condidence intervals
ymax = estimate+(SE*1.96), width = 0.1))+
geom_point(size = 4, shape =21)+
scale_fill_manual(values = c("magenta", "#80ED99"))+
labs(x = "Contrast estimate", y = "Contrast")+
theme(legend.position = "none")+
annotate(x = 1, y = -0.3, label = "***", geom = "text")
contrasts
In the rendered markdown, the labels overlap. This is not the case in the saved image.
ggarrange(estimates, contrasts)
ggsave("figures/modelplots.jpg", height = 4, width = 9, dpi = 400)
3D plots can help visualise data with multiple continuous dimensions. The data set with movement synchrony, movement complexity and dyad-level measures of body competence is perfect for this example.
head(poses_traits1)
## # A tibble: 6 × 6
## # Rowwise:
## ID M_similarity M_entropy mean_BCQ BCQ_half Role
## <int> <dbl> <dbl> <dbl> <fct> <chr>
## 1 114 0.888 0.0555 9.5 High Following
## 2 114 0.907 0.149 9.5 High Following
## 3 114 0.629 0.107 9.5 High Following
## 4 114 0.841 0.0384 9.5 High Following
## 5 114 0.735 0.0637 9.5 High Following
## 6 114 0.557 0.0874 9.5 High Following
The continuous color scale for one axis can help give the plot depth, and thereby help interpretation when rotating it. Click on the plot and drag it around to see different orientations of the axes.
poses_traits1 %>%
plot_ly(x = ~M_similarity, y = ~mean_BCQ, z = ~M_entropy,
marker = list(color = ~M_similarity,
colorscale = c('#FFE1A1', '#683531'),
showscale = TRUE)) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'Synchrony'),
yaxis = list(title = 'Dyadic Body Competence'),
zaxis = list(title = 'Complexity')),
annotations = list(
x = 1.13,
y = 1.05,
text = 'Synchrony',
xref = 'paper',
yref = 'paper',
showarrow = FALSE
))
Alternatively, by creating a new column with which to divide the data into categories, we can help out interpretation too.Click on the plot and drag it around to see different orientations of the axes.
poses_traits1 %>%
plot_ly(x = ~M_similarity, y = ~mean_BCQ, z = ~M_entropy,
color = ~BCQ_half, colors = c("#89398a", "#398a89"), opacity = 0.5) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'Synchrony'),
yaxis = list(title = 'Dyadic Body Compentence'),
zaxis = list(title = 'Complexity')))